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ABSTRACT 

We present a determination of the luminosity functions of massive young stel- 
lar objects (MYSOs) and compact (C) Hll regions within the Milky Way Galaxy 
using the large, well-selected sample of these sources identified by the Red MSX 
Source (RMS) survey. The MYSO luminosity function decreases monotonically 
such that there are few with L >1O^L0, whilst the CHll regions are detected 
up to ~1O^L0. The lifetimes of these phases are also calculated as a function of 
luminosity by comparison with the luminosity function for local main-sequence 
OB stars. These indicate that the MYSO phase has a duration ranging from 
4x10^ yrs for lO^L© to ~7xl0'^ yrs at IO^Lq, whilst the CHll region phase 
lasts of order 3x10^ yrs or ~3-10% of the exciting star's main-sequence life- 
time. MYSOs between IC^Lq and ~1O^L0 are massive but do not display the 
radio continuum or near-IR H I recombination line emission indicative of an H ll 
region, consistent with being swollen due to high ongoing or recent accretion 
rates. Above ~1O^L0 the MYSO phase lifetime becomes comparable to the 
main-sequence Kelvin-Helmholtz timescale, at which point the central star can 
rapidly contract onto the main-sequence even if still accreting, and ionise a CH ll 
region, thus explaining why few highly luminous MYSOs are observed. 
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Introduction 



Massive stars (M > SM©) dominate feedback processes in galaxies due to their strong 
ionising radiation, outflows, stellar winds and supernovae. They have the potential to both 
disrupt their natal molecular cloud and/or trigger further star formation around them. Re- 
gions of massive star formation are usually more heavily embedded within their parent 
molecular cloud and rarer than sites where only low mass stars form, thus we have to look 
further to find statistically sig nificant samples . In addition, massive stars form and evolve 
faster than low mass stars (e.g. IShu et al.lll987l ). These issues mean that there are few large, 
well-selected samples of young massive stars, making it difficult to see global trends over the 
variations of individual sources, and thus restricting our understanding of the formation and 
early evolution of massive stars. Original catalogues of mid-infrared bright but radio-quiet 
massive young stellar objects (MYSOs), where ma.ior accretion is probably still ongoing 

3ut an H II region h as yet to form. 



(see e.g. 



Palla fc Stahler 1993; Behrend fc Maedei 



contained only of order 50 sources f lWvnn-Wilhamslll982[ 



2001} 



pies based on IRAS colours have been constructed ( iPalla et al 



Henning et al. 



1991 



1984). MYSO sam- 



Molinari et al. 1996 



Sridharan et al.ll2002l ) but are often hampered by confusion in regions of high source density 
due to the IRAS satellite's low spatial resolution. Large Hll regions were not included in 
those lists based on comparison with single dish radio surveys, but more compact H ll regions 
were found to contaminate these samples when followed up with higher resolution radio ob- 
servations (Molinari et al. 1998). These criteria therefore tended to bias these samples away 
from complex regions where many massive stars form. 

In order to alleviate this bias, we have conducted a survey to find MYSOs in the Galaxy. 
The Red MSX S ource (RMS) su r vey i s based on a colour-selected sample of ~2000 candi- 
date MYSOs bv iLumsden et al.l (120021) from the 18" resolution MSX mid-IR point source 



catalogue (PSC, lEgan et al. 



1999 



20031 ). The MSX PSC has a peak and mean source den- 



sity of 1812 sources/deg^ and 25 sources/deg^ respectively, so source confusion is not a 
significant issue in the original sample selection. A series of complementary follow-up obser- 
vations have been carried out designed to identify contaminants (low mass YSOs, evolved 
stars, proto-planetary nebulae and planetary nebulae) and provide the characteristics of the 
MYSOs, as well as the ultra-compact and compact Hll regions also present in the sample 
( jUrquhart et al.ll2008b( ). These have included ~1" resolution radio continuum observations 
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( lUrquhart et al.ll2007al . l2009l ). ~1" resolution mid- IR observations (e . g.lMottram et al.ll2007l) 



line observations (jUrquhart et al.ll2007bl l2008al. l201l[) to obtain k inematic distances and low 



outside the SPIT ZER GLIMPSE survev regio ns Jchurchwell et ali boogl). ^=^C0 molecular 



spectral-resolution near-IR spectroscopy (e.g. IClarke et al.ll2006[ ). The spatial resolution of 
the sample is therefore not that of MSX (18") but of the follow-up observations (1—2"). 
Source multiplici ty within the MSX beam is easily revealed by the higher-resolution mid-IR 
follow-up (~25% iMottram et"aDl2007[ boilh . 



Each source within the RMS sample has been classified into one of several categories in 
the database at www.ast.leeds.ac.uk/RMS. 'MYSO's are radio-quiet mid-IR point sources 
which are associated with ^^CO emission. 'Hll regions' are radio-loud and/or resolved in 
the mid-IR indicating emission from Lya heated dust. The division betwee n radio-loud H ll 
region s and radio-quiet MYSOs is clear from the plot of Lrad/Lboi vs. Lboi in lHoare fc Franco 
( 120071 ) . Given the requirement to be an MSX poi nt source (<18"), our H ll regions fall into 
the category of compact and ultra-compact (e.g. iHabing fc Israeli 119791 ) . although we will 
refer to them all as compact from now on for brevity. Compact H ll regions excited by stars 
with spectral types later than Bl may have radio surface brightnesses that are too low to 
be detected in our 1" radio o bservations, but would still exhibit strong hydrogen recombi- 
nation emission with Case-B ( iBaker fc Menzellll938l ) line ratios in our near-IR spectroscopic 
observations (see 



In iMottram et al.l (1201 if ) we used far-IR fluxes from iMottram et al.l (l2010f) and other 



data, along with the spectral energy distribution (SED) model fitter of iRobitaille et al. 
( 120071 ) to obtain luminosity measurement s for ~1200 of these s ources. In this paper we use 



the luminosity distributions obtained by IMottram et al.l (120111 ) to calculate the luminosity 
functions for MYSOs and CH ll regions in the Galactic plane of our Galaxy. These calcula- 
tions will be discussed in §21 the results of which will be presented in §31 We will then use 
these results to obtain estimates for the lifetime of these pha ses as a function of luminosity 
using a method similar to that of IWood fc Churchwelll (Il989[ ). which will be discussed in ^ 
Our conclusions are then presented in ||5l 



2. The Luminosity Function 

We wish to derive a luminosity function defined as the number of objects per unit volume 
per unit luminosity interval in the Galaxy from the observed luminosity distributions, shown 
rebinned to 0.2 dex in Figure 1(a) The first task is to examine the completeness of the sample 
as a function of luminosity. We do this by calculating the volume of the Galaxy within which 
MYSOs and CH lis of luminosity L are detected by the RMS survey. A 3-D grid of cells was 
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set up in cylindrical polar coordinates (i.e. r, 9 and z), centred on the Galactic centre. The 
volume of each cell was calculated {r5r5z59), along with the Galactic coor dinates {^ &: h) 
using the distance between the Earth and the Galactic centre of Rq = 8.5 kpc (IMarshall et al. 



20061 ). The distance between each cell and the Earth was also calculated, which was then 
used to obtain the total flux (Fboi) that would be observed for a hypothetical source of 
luminosity L within the cell. 

Since the RMS sample is effectively a flux- limited sample in the MSX 21 fim band we 
need to convert Fboi to F21 for the hypothetical sources using the observed colours in the RMS 
sample. The obs erved distribution of F hni / -F21 was obtained from SED model fits for young 



RMS sources by iMottram et al.l (120 111 ), and takes the for m of a Gaussian in l og(Fhoi/-F2i) 
with mean 1.43 and standard deviatio n 0.27 (see fig 7. , iMottram et al.l 120111). It is wel l 



known the shape of the SEDs of CH lis fIChini et al.lll987n and MYSO flHenning et al.lll990l ) 
are on average very similar, and no significant variation in -Fboi/-F2i between MYSOs and 
Hll regions was found in the model SEDs. The mean ratio was therefore used to calculate 
the MSX 21 fim fiux in Wm~^ from Fboi for a hypothetical source of luminosity L in the 
cell. These are then divided by the bandwidth of the MSX 21 /im band in order to convert 
the fiux to Janskys. 

The completeness of the MSX survey at 21 fim varies with i and b, particularly between 
the inner and outer galaxy due to the number of times each region of the sky was scanned, 
so we implement different flux limits (-Fumit) for these two regions. Given that most RMS 
sources lie in the Galactic mid-plane, we use F^^it = 3.5 Jy for 270° > i > 90° and 4.25 Jy 
for 90° > i > 270° based on the 50% completion limits calculated for | & | < 1° by Davies et 
al. (2010, in prep). If a source with luminosity L within the cell has F21 > Fn^a then the 
volume of the cell is added to the total volume included in the RMS survey. 

The observed distributions of young RMS sources with luminosities L > 8.2x10'^Lq 
(discussed in Q are shown in Figures 1(b) and 1(c) in terms of galactocentric radius (r) 
and distance from the Galactic mid-plane (z) respectively. Given that these are far from 
uniformly distributed within the Galaxy, a weighting function was used to weight the con- 
tribution made to the total volume by each cell according to the observed source distribu- 
tion. This has the functional form of a thin disk with a hole in the centre as modelled by 
Robin et al.l (120031 ) for the Galactic stellar population of the thin disk. We adopted length 
scales for the disk (rd) and central hole (rh) consistent with our observed distribution in 
galactocentric radius (see Figure 1(b)) and also allowed for an offset (rg 



from the Galac- 



tic centre to account for the observed lack of sources with r < 2 kpc. The distribution of 



sources in terms of z (Figure 1(c) ) is that of an exponential with scale-height z^^. The weight 
is therefore given by: 
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Fig. 1. — Luminosity function inputs, (a): the number of s ource s vs. luminosity for RMS 
MYSOs (black) and CHll regions (red) from lMottram et al.l (l201ll ) rebinned to 0.2 dex. The 
black arrow indicate the luminosity of a BOV star from Table [H (b,c): the distribution of 
luminous (L > 8.2x10^Lq) sources in the RMS survey as a function of galactocentric radius 
(r) and distance from the Galactic mid-plane {z) respectively. The red lines indicate the 
least-squares fit of the weighting function (equation [1]) to the data. 
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weight{r, z) 



exp 



— z 



exp 



exp 



-{r - To) 



A least-squares minimisation was used to obtain the values of the parameters in this function 
from the distribution of RMS sources, which resulted in = 6.99 ± 0.22 kpc, r-^ = 1.71 ± 0.23 kpc, 
ro = 2.20 ± 0.10 kpc and z^ = 39.0 ± 1.8 pc, shown by the red solid lines in Figures 1(b) 



and 1(c) , with x^g^=3.8 and 1.2 for the r and z fits respectively. This result does not change 
within the errors if luminosity cuts corresponding to 10% or 90% completeness are used 
instead of 50%. The asymmetry in z and the slight increase in the number of sources at 
r = 10 kpc are prob ably due to local features near the Sun (i?o = 8.5 kpc, Zq = 15 pc, 
Marshall et al.ll2006l ). Both the shape of the distributions and the fit results ar e simil ar to 
those obtained pre viously for site s of na assive star formation by iBronfman et al.l ( 120001 ) and 
for Hll regions by iPaladini et al.l (120041 ). 



The total weighted volume for which a source with luminosity L should be included 
within the RMS survey is therefore calculated by integrating th e weighted cell volurn e over 



all cells with F21 > Fn^it, -5° | 6 | < 5° and 10° > £ > 350° (see lLumsden et al.ll2002h . This 
process was then repeated for logarithmic steps in luminosity. Increments of 0.1 kpc, 0.5° 
and 0.01 kpc were used in r, 6 and z respectively, and all values were calculated for the cell 
centre. 

The error in the luminosity is taken to be the bin width (0.2 dex) and the error in 
the luminosity function is calculated from statistical uncertainty using Poissonian statistics. 
Only bins with at least 5 sources in them are considered in the follo wing analysis, ensuriri g 
that while the average uncertainty in source luminosity is ~34% ( iMottram et al.l 1201 ll ). 
dominated by our conservative estimate of distance error, the ensemble uncertainty is smaller 
than the bin size. 



3. Results 



The volume of the Galaxy sampled by the RMS survey at each luminosity bin is shown 
This indicates that the RMS survey should be > 50% complete above 



in Figure 2 (a 



~8.2x 1O^L0, which corresponds to a spectral type of ~B0.5V using the scale presented in 
Table [TJ We restrict the source distribution calculations to above the 50% completion limit, 
as below this luminosity the volume correction required is > 2. The luminosity functions were 



then obtained by dividing the source numbers in Figure 1(a) by the volume in Figure 2(a^ 
and converting to per dex in L. 
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Fig. 2. — Luminosity function results, (a): Volume of the Galaxy within which sources should 
be detected in the RMS survey vs. luminosity. The green and blue dashed lines indicate 
100% (13.6 kpc^) and 50% and completion respectively, (b): The luminosity functions for 
MYSOs (black) and CHll regions (red), (c): The combined luminosity functions for sources 
identified as either MYSOs or CHll regions in the RMS survey. The results of the total 
least-squares fits to the luminosity functions over limited ranges of L are indicated by black 
and red solid lines, and are discussed in ^ The spectral type scale is presented in Table [1] 
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The resulting luminosity functions are shown above the 50% completion limit in Fig- 
ure 2(b) Over the luminosity range 10^ < L < 10^ Lq, the MYSO luminosity function 
falls steadily such that log(<^MYSo) = (-1-3 ± 0.2) log(L/L0) + (6.9 ± 0.9). Consequently, 
no significant population of MYSOs is observed above 10^ Lq or earlier than ~06.5V. 

There is a flattening in the CHll region luminosity function around L = 4x10^ Lq 
(09V). Above this the luminosity function is fitted by log(0cHii) = (-0.8 ± 0.1) log(L/LQ) + 
(5.0 ± 0.8). Although the turnover in slope of the CHll region luminosity function is below 
the 100% completeness limit, the fact that there is no similar break in the MYSO luminosity 
function that continues to lower luminosities means that the break is probably real. 



Figure 2(c) shows the luminosity function for the combined sample of all MYSOs and 
CHll regions identified in the RMS survey. This shows the overall function for the young 
mid-IR bright phase of massive star formation regardless of whether the source is ionising 
its surroundings or not. Above L = 4x10^ Lq this combined luminosity function is fitted by 

log(0Combined) = ("0.9 ± 0.2) log(L/L0) + (5.8 ± 0.8). 



4. Timescales 



Wood fc Churchwelll (Il989l ) used the number of main-sequence stars to estimate the 



lifetime for all UGH ll regions, independent of luminosity, from the number of such sources 
they observed. The same method can be used to estimate the lifetimes of MYSOs and CH ll 
regions as a function of luminosity, using (j)MYSo{L) and (j)cmi{L), i.e.: 



^MYSo(-^) — ^Ms(-^) 



4>Ms{L) 



(2) 



where the main-sequence lifetime (^ms) of an OB star is given by log(tMs/yrs) = (11.4 ± 1.2) - 
(1.5 ± 0.5) log(L/L0) + (0.11 ± 0.05) log^fL/L^) calculated from a l east-s quares fit to the 
results of the models including stellar rotation by iMeynet fc Maeded ( 120031 ). These models 
predict the correct ratio of red to blue supergiants, and are the only models which extend 
beyond 60 Mq. For M>2OM0 we use the O star lifetime, while for M<2OM0 we use the 
hydrogen-burning lifetime. 

The local main-sequence OB star lu minos i ty fun ction in terms of My was obtained from 
re-binning the results on dwarf stars of iReedl (120051 ) corrected to the My scale in Table (H 



resulting in a least-squares fit of the form log(0oB(^v) /stars kpc ^mag ^) = (4.12 ± 0.13) 
+ (0.60 ± 0.06) My. This is shown in Figured 



- 9 - 




Fig. 3. — Dwarf OB star luminosity function, calculated from the results of iReedl (120051 ) 
using the My-Spectral type scale given in Table [1] 
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This was then converted from (/)oB(Afv) to (j)oB{L) using My = (1.5 ± 0.3) - (0.6 ± 0.1) log(L/L0) 
(0.10 ± 0.02) log^(L/L0) which is a least-squares fit to the data in Table [H We mul- 
tiply our luminosity functions by the weight for the Galactic location of the Sun (i.e. 
weight (i?0 5-^o) = 0.302 from Eqn. [T]) before using Eqn. [2] so that they can be compared 
directly to the local OB luminosity function. 

The MYSO and H ll region lifetimes obtained as functions of luminosity are shown in 
Figure m The difference between the fits to the data for My vs. log(L), 0ob vs. My and 
log(tMs) vs. log(L) and the data themselves on the timescale calculations result in a system- 
atic uncertainty of approximately a factor of 10, as shown by the black error bar. Fits to the 
data give log(tMYSo/yrs) = (-0.80 ± 0.24) log(L/LQ) + (8.9 ± 1.1) and log(tcHii/yrs) = (- 
0.13 ± 0.16) log(L/L0) + (6.1 ± 0.8). The systematic error in the slope is of order 0.5, 
obtained by considering the errors in each step of the transformation. 

The timescale of MYSOs decreases by a little less than an order of magnitude be- 
ween lO'^L© and lO^L^r^, while the CHl l region lifetime is flat with respect to luminosity at 



~3xl0^yrs. IWood fc Churchwelll (119891 ) concluded that the lifetimes of Hll regions powered 
by O stars are about 15% of their main-sequence lifetimes. Here we find the CHll region 
lifetimes are ~3— 10% of their main-sequence lifetimes, a lthough the main-seq u ence l ifetimes 



we are using are about twice as long as those used by IWood fc Churchwelll (119891 ) so our 
results are broadly consistent with theirs. 



5. Discussion & Conclusions 

We have obtained the luminosity functions of MYSOs and compact H ll regions in the 
Galactic plane. These clearly show that there is no significant population of the mid-IR 
bright, radio-quiet MYSOs above ^lO^L©. This lack of high luminosity MYSOs could mean 
that the immediate progenitors of the most luminous CH ll regions are mid-IR faint or dark. 
Luminous far-IR cores and methanol masers that are too faint to be in our MSX-selected 
sample are being found (e.g. 



Chambers et al.ll2009l : lEUingsenl 120061 ) 



However, it is also important to consider the luminosity function in terms of the mech- 
anism that must inhibit the formation of H ll regions around MYSOs, which are energetic 
enough to ionise their surroundings if they were on the main-sequence. Previously it has 



Walmslev 


1995; 


Keto 


2003 



Tan Sz McKedl2003l ). In this case the H ll region would have a high enough emission measure 
(EM) to be optically thick in the radio continuum, e.g. Ti cm = 3x10^^° EM (cm"^pc)^l 
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Tabl e 1: Adopted propert ies o f main-sequence stars. The surface temperatures are taken 
(120051 ) and iBoehm-Vitensd ( 1l98ll ) for O and B stars respectively. The 



from Martins et al. 



luminosities and masses were calculated using a 1 Myr isochrone, chosen to better represent 



th e observed local OB pop ulation, from the models with stellar rotation 



of iMevnet fc Maeded (12003 



'iiu; 



=300 kms" 



My was calculated using the equations given in iMartins et al. 



( 120051 ) and iLanz fc Hubenyl (120071 ) for O and B type stars respectively. The main-sequence 
lifetimes were calculated from the quadratic fit (see Q to the main-sequence lifetimes of 



Mevnet & Maeder 


(2003) 


see 


Martins et al. 


2OO5I). 
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Fig. 4. — The lifetime of MYSOs (black) and CH ll regions (red) vs. luminosity. The red 
and black dashed lines indicate the results of least-squares fits to the functions over limited 
ranges of L and are discussed in ^ The green and blue solid lines indicate the main-sequence 
Kelvin- Helmholtz timescale and main-sequence lifetime respectively. The black vertical error 
bar in the top left corner indicates the systematic uncertainties. 
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( 10sterbrocklll989l ) . Unless such regions have EM > 10 near-IR recombination hnes will still 



be optically thin since the ir line centre optical depths are much lower, e.g. t^^-^ ~2 x 10 ^ ti 



( iHummer fc Storeyl 119871 ). This would lead to narrow lines with high equivalent widths 
whereas in MYSOs where we detect the NIR continuum, we actually observe broad lines 
with low equivale nt widths consistent with optically thick emission from stellar winds (e.g. 
Bunn et al.lflQQsh . 



Alternatively, as first suggested by lHoare fc Francol (120071 ). the MYSO could be too cool 
to ionise its surroundings on account of being swollen due to high accretion rates. Recent 



stellar evolution calculatio i is including high accr etion rates ( lYorke &: Bodenheimerl 12008 



Hosokawa fc Omukail l2009l : iHosokawa et al.l |2010[ ) have suggested that the radius of stars 
with masses between ~8M0 and ~3OM0 increase dramatically (i.e. of order 10-100 times 
the main-sequence size) due to high accretion rates. Formation of an Hll region around 
an actively accreting star only begins when it is massive enough that the Kelvin-Helmholtz 
timescale is less than the accretion timescale. We have shown that the timescale associated 
with MYSOs becomes comparable with the main-sequence Kelvin Helmholtz timescale at 
~1O^L0. That we see no significant population of MYSOs with L >1O^L0 is consistent with 
this idea. For the most massive MYSOs, i.e. those that have reached ~3OM0, the accretion 
rate must be between 1O~^M0 yr~^ and 1O~^M0 yr~^, while the average accretion rates for 
less massive MYSOs are likely to be lower given the lifetimes we have obtained. 

The luminosity functions of MYSOs and H ll regions presented in this paper can be used 
to test models of the time-evolution of the mass accretion rate and stellar properties. Such 
modelling will be discussed in a future publication (Davies et al. 2010, in prep.). 
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